# File: BL.R
# R source code for computing Black-Litterman Model
# Author: Congxing Cai (congxing@stanford.edu)
black.litterman <- function(data, P, Q) {
	# Compute the posterior return using black-litterman model
	
	# Step 1: Compute equilibrium risk prima
	E <- get.covariance(data)
	cap.wt <- get.market.capital.weight(data)
	d <- compute.risk.aversion()
	II <- reverse.optimization(d, E, t(as.matrix(cap.wt)))
	II <- t(as.matrix(II)) # convert to N*1 vector
	
	# Step 2: Compute posterior mu and sigma
	tau <- 0.02

  	# omega for absolute view matrix, P
  	omega <- diag(diag(E))*tau # Based on the computation described in the report
	mu.bl <- solve(solve(tau*E) + t(P) %*% solve(omega) %*% P) %*% (solve(tau*E) %*% II + t(P) %*% solve(omega) %*% Q)
	sigma.bl <- E + solve(solve(tau*E) + t(P) %*% solve(omega) %*% P)

  	return(list(mu=mu.bl, sigma=sigma.bl, aversion=d))  
}

reverse.optimization <- function (d, E, w) {
	# Reverse optimization: II = d * E * w
	# Args:
	#	 d  - risk aversion parameter
	#  	 E  - covariance matrix n x n based on individual assets in portfolio
	#  	 w  - weights of portfolio assets based on market capitalization, n x 1 vector
	# Returns:
	#   II - implied returns based on n asset in a portfolio	
	II <- data.frame(t(d * as.matrix(E) %*% w))
  	return(II)
}

compute.risk.aversion <- function(risk.free.rate = 0) {
  # The implied risk aversion coefficient can be estimated by dividing the expected excess return by the variance of the portfolio. We use S&P500 as the benchmark.
  m_SP500<-read.table("../data/m_SP500.csv", header=TRUE, sep=",")
  # this returns are log returns, so must be aware to compare vs log returns
  benchmark<-log((m_SP500$Adj.Close[2:length(m_SP500$Adj.Close)] / m_SP500$Adj.Close[1:(length(m_SP500$Adj.Close)-1)]))
  lambda = mean((benchmark) - (risk.free.rate)) / var((benchmark))
  return( as.double(lambda) )
} 

get.market.capital.weight <- function(data) {
	# Compute the market capitalization weights w.
	# Args:
	#   data: the data frame containing the msize and nfirms for smlo, smme, smhi, bilo, bime, bihi 
	# Returns:
	#	FFw market capitalization weights.
	data <- data[nrow(data),]
	msize <- data.frame(data$date, data$smlo_msize, data$smme_msize, data$smhi_msize, data$bilo_msize,	data$bime_msize, data$bihi_msize)
  	nfirms <- data.frame(data$date, data$smlo_nfirms, data$smme_nfirms, data$smhi_nfirms, data$bilo_nfirms, data$bime_nfirms, data$bihi_nfirms)
  	cap	<- data.frame(date=data$date, msize[,-1] * nfirms[,-1])
  	cap <- data.frame(cap, total.mkt.cap = rowSums(cap[,-1]))
 	 # dividing each FF portfolio's mkt cap by sum of market caps.
  	cap.weight <- data.frame(cap[,2:7] / cap$total.mkt.cap)
  	names(cap.weight) <- c("smlo.mkt.wt","smme.mkt.wt","smhi.mkt.wt","bilo.mkt.wt","bime.mkt.wt","bihi.mkt.wt")
  	return(cap.weight)
}

get.covariance <- function(data) {
	# Compute the covariance matrix based on individual assets in portfolio.
	# Args:
	#	data: the data frame containing the vw_ret for smlo, smme, smhi, bilo, bime, bihi
	# Returns: 
	#   The covariance matrix.
	ret <- data.frame(data$smlo_vwret, data$smme_vwret, data$smhi_vwret, data$bilo_vwret, data$bime_vwret, data$bihi_vwret)
	names(ret) <- c("smlo", "smme", "smhi", "bilo", "bime", "bihi")
  	E <- cov(ret)
  	return(E)
}